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Abstract 

EUPDF is an Eulerian-based Monte Carlo PDF solver developed for 
application with sprays, combustion, parallel computing and unstructured 
grids. It is designed to be massively parallel and could easily be coupled 
with any existing gas-phase flow and spray solvers. The solver accommodates 
the use of an unstructured mesh with mixed elements of either triangular, 
quadrilateral, and/or tetrahedral type. The manual provides the user with 
the coding required to couple the PDF code to any given flow code and a basic 
understanding of the EUPDF code structure as well as the models involved 
in the PDF formulation. The source code of EUPDF will be available with 
the release of the National Combustion Code (NCC) as a complete package. 
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I. Nomenclature 


A pre-exponential coefficient in an Arrhenius reaction-rate term 
a n outward area normal vector of the nth surface, m 2 
C p specific heat, J/(Kg K) 

c n convection/diffusion coefficient of the nth face, kg/s 
D turbulent diffusion coefficient, m 2 /s 

E a activation energy in an Arrhenius reaction-rate term 
h specific enthalpy, J/Kg 

J-* diffusive mass flux vector, Kg /ms 

k turbulence kinetic energy, m 2 /s 2 

molecular weight of ith species, kg/kg-mole 
droplet vaporization rate, Kg/s 
lk y eff effective latent heat of evaporation, J/Kg 

N av number of time steps employed in the PDF time-averaging scheme 

Nf number of surfaces contained in a given computational cell 

N m total number of Monte Carlo particles per grid cell 

N p total number of computational cells 

rik number of droplets in kth group 

P pressure, N/m 2 

P r Prandtl number 

p joint scalar PDF 

R u gas constant, J/(Kg K) 

r radial coordinate (gas-phase equations), m 

S a liquid source contribution of the a variable 

S c Schmidt number 

S m ic liquid source contribution of the gas-phase continuity equation 

Smie liquid source contribution of the gas-phase energy equation 

Smim liquid source contribution of the gas-phase momentum equations 
S m ts liquid source contribution of the gas-phase species equations 

T temperature, K 

t time, s 

U{ velocity component in the ith direction, m/s 

w a chemical reaction rate, 1/s 

Cartesian coordinate in the ith direction, m 
Yj mass fraction of jth species 

x spatial vector 

greek symbols 
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turbulent diffusion coefficient, kg/ms 
local time step used in the PDF computations, s 
local time step in the flow solver, s 
computational cell volume, m 3 
Dirac-delta function 

species mass fraction at the droplet surface 
rate of turbulence dissipation, m 2 /s 3 
thermal conductivity, J/(ms K) 
dynamic viscosity, kg/ms 
density, kg/m 3 
mole fraction 
dimensionality of ^-space 
independent composition space 


u> turbulence frequency, 1/s 
r stress tensor 


superscripts 

Favre averaging 

time averaging or average based on the Monte Carlo 
particles present in a given cell 
It fluctuations 

subscripts 

a scalar variable of the joint PDF 
f represents conditions associated with fuel 
g global or gas-phase 

i coordinate or specie indices 

j specie indices 

k droplet group or liquid phase 

l liquid phase or laminar 

m conditions associated with N m 

n nth-face of the computational cell 

o initial conditions or oxidizer 

p grid cell 

s represents conditions at the droplet 
surface or adjacent computational cell 
t conditions associated with time 

, partial differentiation with respect 

to the variable followed by it 
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II. Introduction 


The gas-turbine combustor flows are often characterized by a complex 
interaction between various physical processes associated with the interaction 
between the liquid and gas phases, droplet vaporization, turbulent mixing, 
heat release associated with chemical kinetics, radiative heat transfer as- 
sociated with highly absorbing and radiating species, among others. The 
rate controlling processes often interact with each other at various disparate 
time and length scales. In particular, turbulence plays an important role 
in determining the rates of mass and heat transfer, chemical reactions, and 
liquid phase evaporation in many practical combustion devices. Most of the 
turbulence closure models for reactive flows have difficulty in treating non- 
linear reaction rates. 1-2 The use of assumed shape PDF methods was found 
to provide reasonable predictions for pattern factors and NOx emissions at 
the combustor exit. However, their extension to multi-scalar chemistry be- 
comes quite intractable. The solution procedure based on the modeled joint 
composition PDF transport equation has an advantage in that it treats the 
nonlinear reaction rates without any approximation. This approach holds 
the promise of modeling various important combustion phenomena relevant 
to practical combustion devices such as flame extinction and blow-off limits, 
and unburnt hydrocarbons (UHC), CO, and NOx' predictions. 2 

However, in the PDF transport equation, all of the composition variables 
are present as independent variables in addition to the space and time vari- 
ables. Because of the large dimensionality of a joint scalar PDF transport 
equation, a deterministic solution becomes impractical. 2 However, Monte 
Carlo methods are better suited over other numerical methods because of 
the advantage that the computational effort rises only linearly with the di- 
mensionality of the PDF. But they tend to be computationally very time 
consuming and require a large computer memory for application to 3D flows. 
However, the success of any numerical methodology used in the study of 
practical combustion flows depends not only on the modeling and numeri- 
cal accuracy considerations; but its applicability would be dictated mainly 
by the available computer memory and turnaround times afforded by the 
present-day computers. 

Other than some simple cases, there is a limited experience with the 
application of the PDF method based on the solution of the PDF transport 
equation to the calculations involving three-dimensional practical combus- 
tion flows. 3 With the aim of demonstrating its viability to the modeling of 
practical combustion flows, we have undertaken the task of extending this 
technique to: (1) the modeling of sprays in order to facilitate simulation of 
gas-turbine combustion flows, (2) parallel computing in order to facilitate 
large-scale combustor computations, and (3) unstructured grids in order to 
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facilitate computations on complex combustion geometries. Also, several nu- 
merical techniques are outlined for overcoming some of the high computer 
time and storage limitations associated with the Monte Carlo simulation of 
practical combustor flows. 

Some of our previous work on the Monte Carlo simulation could be found 
in Refs. 4 to 7. Initially, the emphasis of our work was on extending the joint 
scalar Monte Carlo PDF method to the modeling of compressible reacting 
flows. 4 The Monte Carlo solver was used in conjunction with several density- 
based codes for the mean-velocity and turbulence fields. Several averaging 
procedures introduced in Refs. 4 and 5 proved to be useful in providing 
smooth Monte Carlo solutions to the CFD solver. The PDF method provided 
favorable results when applied to several supersonic diffusion flames. 4-5 

Later on this approach was further extended to the modeling of spray 
flames and parallel computing. 6 This method combined the novelty of the 
PDF method with the ability to run on parallel architectures. This algo- 
rithm was implemented on the Cray T3D at NASA Lewis Research Center, 
a massively parallel computer, with an aggregate of 64 Processor Elements 
(PEs). The computer code was written in Cray MPP (Massively Parallel 
Processing) Fortran. The application of this method to both open as well as 
confined axisymmetric swirl-stabilized spray flames showed good agreement 
with the measured data. 6 Preliminary estimates indicated that it was well 
within modern parallel computer’s capacity to do a realistic gas-turbine com- 
bustor simulation on grid sizes of 125,000 nodes and a total of 12.5 million 
Monte Carlo particles with reasonable turnaround times. 

It is well known that considerable effort usually goes into generating 
structured grid meshes for gridding up practical combustor geometries which 
tend to be very complex in shape and configuration. The grid generation 
time could be reduced considerably by making use of existing automated 
unstructured grid generators. 8 With the aim of advancing the current multi- 
dimensional computational tools used in the design of advanced technology 
combustors, we have recently extended this technique to unstructured grids 
following the guidelines established for the development of the National Com- 
bustion Code (NCC). 7 NCC is being developed in the form of a collaborative 
effort between NASA LeRC, aircraft engine manufacturers, and several other 
government agencies and contractors. 9 Some of the salient features of our 
work in Ref. 7 are summarized below: 

(1) The scalar Monte Carlo PDF solver has been extended in a cell- 
centered, finite volume context to mixed unstructured grid elements of tri- 
angular, quadrilateral, tetrahedron, wedge, and hexahedron elements. 

(2) The PDF code was rewritten in Fortran 77 with PVM calls for par- 
allel computing which enables the computations to be performed with equal 
ease on both massively parallel computers as well as workstation clusters. 
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(3) The PDF module was coupled with Pratt and Whitney’s 
CORSAIR 10 - an unstructured flow solver, and LSPRAY 11 - a Lagrangian 
spray solver, which were selected to be the integral components of the NCC 
cluster of modules. LSPRAY was developed for application with unstruc- 
tured grids and parallel computing. 

(4) Our experience has shown that by introducing the concept of local 
time-stepping into the Monte Carlo PDF computations, we were able to 
avoid the occurrence of the so-called frozen condition 3 and, thereby, able to 
compute the solution with the use of a relatively fewer number of stochastic 
particles in flows characterized by recirculation, swirl, and boundary layers. 

(5) The effect of temperature variation was taken into account in the 
evaluation of both temperature and specific enthalpy of a gaseous mixture. 

(6) The PDF solver receives the mean-velocity and turbulence fields from 
the flow solver and the source terms arising from the liquid-phase contribu- 
tion from the spray solver. It in turn provides the species and temperature 
solution to the spray and CFD solvers. 

The PDF method seems to provide favorable agreement when applied 
to several supersonic diffusion flames as well as several other swirl-stabilized 
spray flames. 4-7 


III. Composition Joint PDF Equation 

The transport equation for the density- weighted joint PDF of the com- 
positions, p , is: 


[pp),t + [pUip],n + [pu> a (g)p\,4, a = 

{Mean convection } {Chemical reactions} 

-\p < U " I t > PU ~[p< - I ± > PUa 

{Turbulent convection} {Molecular mixing} 

- [p < -s a I > p\^ a 
p ~ 

{Liquid — phase contribution} 


( 1 ) 


where 

w a = 

< u ■' | 3 /; > = 

<\J°,A±> = 


chemical source term for the o-th composition variable, 
conditional average of Favre velocity fluctuations, 
conditional average of scalar dissipation, and 
conditional average of spray source terms. 
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The terms on the left hand side of the above equation could be evalu- 
ated without any approximation but the terms on the right hand side of the 
equation require modeling. The first term on the right represents transport 
in physical space due to turbulent convection. 2 Since the joint PDF, p , con- 
tains no information on velocity, the conditional expectation of < u” | 0 > 
needs to be modeled. It is modeled based on a gradient-diffusion model with 
information supplied on the turbulent flow field from the flow solver. 2 

- < u" | 0 > p = (2) 

The second term on the right hand side represents transport in the 
scalar space due to molecular mixing. A mathematical description of the 
mixing process is rather complicated and the interested reader is refered to 
Ref. 2. Molecular mixing is accounted for by making use of the relaxation 
to the ensemble mean submodel 1 as it seemed to provide rather satisfactory 
results. 12 


< -J?* I </> >= -C>(0 a - 0 a ) (3) 

P 

The third term on the right hand side represents the contribution from 
the the spray source terms. The conditional average is modeled based on the 
average values of species and enthalpy: 

< -s a I 0 >= — TT7 y2n k m k (e as - <j> Q ) (4) 

p — pAV 

where <j) a = Y Qy a = 1 , 2 , ..., s — a — \ 

< -s Q I 0 >= — f— y n k m k {-l kteJ} + h ks - <f> a ) (5) 

p — p A V 

where <ft a = h. 

where e QS is a mass fraction of the evaporating species at the droplet surface. 
Here we assumed that the spray source terms could be evaluated independent 
of the fluctuations in the gas phase compositions of species and enthalpy. Eqs. 
4 to 5 represent the modeled representation for the conditional averages of 
the spray contribution to the PDF transport equation. 

IV. Solution Algorithm 

In order to facilitate the integration of the Monte Carlo PDF method 
in a finite-volume context, the volume integrals of convection and diffusion 
in Eq. (1) were first recast into surface integrals by means of a Gauss’s 
theorem. 3 Partial integration of the PDF transport equation would yield: 
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p p (±, t + At) = (l- j^)p P {±, t) + J2 0 

- A t[w Q (f)p}^ a - Af[< I V> > p} yi>a - A t[< y a \±> p],v, Q (6) 

with subscript n refers to the nth-face of the computational cell. The coeffi- 
cient c n represents the transport by convection and diffusion through the nth 
face of the computational cell, p. The convection/diffusion coefficients in the 
above equation are determined by one of the following two expressions: 

° n = ^ A^f Af/J + max l°’ -P^n-U-n] 

2(2 CL 

c n = max[\0Apa n .u n \,T4 -^=^^ ) ] - 0.5pa n .u„ 

and 

Cp — ^ ' c n 

n 

In both the above expressions for c n , a cell-centered finite- volume derivative 
is used to describe the viscous fluxes; but an upwind differencing scheme is 
used for the convective fluxes in the first expression and a hybrid differencing 
scheme in the second. 

Numerical Method Based on Approximate Factorization 


The transport equation is solved by making use of an approximate fac- 
torization scheme. 2 Eq. 6 can be recast as 


p p (rp,t + At) = 

(/ + A tR)((I + AtS)(I + A tM)(I + AtT)p p {±,t) + 0{At 2 ) (7) 

where I represents the unity operator and T, M, S, and R denote the opera- 
tors associated with spatial transport, molecular mixing, spray, and chemical 
reactions, respectively. The operator is further split into a sequence of inter- 
mediate steps: 


Pp(tM) = (/ + AtT)p p (^,t) 


(8) 
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P"(±,t) = (i + AtM)p;(±,t) 


( 9 ) 


g“(0,<) = (/ + AiS)^*(0,f) (10) 

p r (l,t + At) = (I + AtR)p^(±,t) (11) 

The operator splitting method provides the solution for the transport of p 
by making use of a Monte Carlo technique. In the Monte Carlo simulation 
the dentity weighted PDF at each grid cell is represented by an ensemble 
of N m stochastic elements where the ensemble-averaged PDF over N m delta 
functions replaces the average based on a continuous PDF. 2 

=< p P w >= 4- e w - n < 12 > 

ly ™ n=l 

The discrete PDF p pm (t, [>) is defined in terms of N m sample values of 4> n , 
n = 1,2,3 ...N m . The statistical error in this approximation is proportional 

toA, ” /2 - 

Using the operator splitting method, the solution for the PDF transport 
equation is obtained sequentially according to the intermediate steps given 
by Eqs. 8 to 11. 

Convection/Diffusion Step 


The first step associated with convection/diffusion is given by: 


P* p {±, <) = {I + A tT)p p (rp, t) = 

+ ( 13 ) 

This step is simulated by replacing a number of particles ( = the nearest 
integer of Cr A^ m ) at cj)p(t) by randomly selected particles at (f> n (t). 

Numerical Issues Associated With Fixed Versus Variable Time Step 
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It is obvious from the above equation that a necessary criterion for 
stability requires satisfaction of < 1. This restriction needs to be 

satisfied for all grid nodes at the same time. However, this criterion tends 
to be too restrictive for applications involving complex combustor chamber 
flows where resolution considerations require concentration of the grid lines 
more in certain regions of the flowfield than the others. For example, more 
grid lines are clustered in regions where boundary layers are formed. The 
time step as determined by a limited region in the computational domain 
can lead to a frozen condition in some other nodes where there may be no 
elements for shifting across the neighboring cells. In order to avoid this frozen 
condition, the following criterion 


c n AtN m 
pAV > 

has to be satisfied at all grid nodes; but such a restriction could invariably 
lead to a prohibitively large cpu time. Scheurlen et al. 3 were the first ones to 
recognize the limitations associated with the use of a fixed time step in the 
Monte Carlo PDF computations. 

However, our experience has shown that this problem can be overcome 
by introducing the concept of local time-stepping which is a convergence 
improvement technique widely used in many of the steady-state CFD com- 
putations. In this approach, the solution is advanced at a variable time step 
for different grid nodes. In our present computations, it is determined based 
on 


At = min(C lf At f , ^ — -) 

* dmlc) 

where Ctj and C\ are constants and were assigned the values of 4 and 2.5, 
respectively, A tf is the local time step obtained from the flow (CORSAIR) 
module, and S m lc = J2 n k^k- The time step is chosen such that it permits 
transfer of enough particles across the boundaries of the neighboring cells 
while ensuring that the time step used in the PDF computations does not 
deviate very much from the time step used in the flow solver. 

Molecular Mixing Step 

The second step associated with molecular mixing is given by 


d<ba 

dt 


4*<x) 


The solution for this equation is updated by: 


( 14 ) 


where u 


C = K + (fa - fa)e~ C *“ A ' 

e/k , and C<f> was assigned a value of 1. 


( 15 ) 
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Spray Step 


The third step associated with the spray contribution is given by 

IT = 

where <f> a — V a ,Q = 1,2 , s = a — 1 


d(pOr 


1 


^ ffc,e// T ^cr) 


( 16 ) 


(17) 


dt pAV 

where ^ = h. 

The solution for the above equations is upgraded by a simple explicit 
scheme: 

2 *** A tY,n k m k AtJ2n k m k 

4> a = rxiT + <M 1 7X77 ) ( lb ) 


pAV 


pAV 


where « < <7 — 1 


<t> 


'k'k'k 

a 


A t^2n k m k 

pAV 


( lk,ef f T ^fcs) 


+ C(1 


At g 

pAV 


(19) 


where a = a. 

After a new value for enthalpy is updated, temperature is determined 
iteratively from the following equation: 


where 


Na-l 

h = y> h ' 

t=i 


( 20 ) 


hi — h°ji + f C p iyidT , 

JT rcf 

C pi = ^(A„- + A 2l T + A 3l T 2 + A 4i T 3 + A 5i T 4 ), 
h°j { is the heat of formation of ith species, and R u is the universal gas constant. 

Reaction Step 


Finally, the fourth step associated with chemical reactions is given by: 
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where <j> a = Yf. 


Wo Af P^S \0.25/' \ 1.25 — ( 4j(M 

<ft " ° p Vy j W 6 

where <f> a = Y 0 . 


(22) 


dp(f>oc 

dt 

where <f> a = h. 

The numerical solution for Eqs 20-23 is integrated by an implicit Eu- 
ler scheme. 13 The resulting non-linear algebraic equations are solved by the 
method of quasi-linearization. 14 



Details of combustion chemistry 


Combustion is modeled by finite-rate kinetics with a single step global 
mechanism of Westbrook and Dryer 15 for n-heptane. This global combustion 
model is reported to provide adequate representation of temperature histories 
in flows not dominated by long ignition delay times. For example, the overall 
reaction representing the oxidation of the n-heptane fuel is given by 

CyH\q + 11(02 + 3.76/V2) — ► 

700 2 + 8 H 2 0 + 41.367V 2 (24) 

Because of the constant Schmidt number assumption made in the PDF 
formulation, based on atomic balance of the constituent species, the mass 
fractions of N 2 , C0 2 , and H 2 0 can be shown to be related to the mass 
fractions of 0 2 and CjH\q by the following expressions: 

Vh 2 o — I< 2 — KiI\ 2 yo 2 — ^2yc 7 H 16 
yco 2 = Tv 2 A 3 — A ! A 2 A ■ i \jo 2 — K 2 I< 3 yc 7 H l6 (25) 

Vn 2 = l- I <2 - A 2 A 3 - yo 2 ( 1 - I<iK 2 - K l K 2 K 3 )- 
vc 7 h 16 ( 1 - A' 2 - AVCj) 

where I<i = 4.29, I< 2 = 0.08943, and I< 3 = 2.138. 

Using Eq. 25 results in considerable savings in computational time as it 
reduces the number of variables in the PDF equation from five (four species 
and one energy) to three (two species and one energy). 

Revolving Time- Weighted Averaging 
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It is noteworthy that although local time stepping seems to overcome 
some of the problems associated with the PDF computations, the application 
of the Monte Carlo method requires the use of a large number of particles 
because the statistical error associated with the Monte Carlo Method is pro- 
portional to the square root of N m which makes the use of the Monte Carlo 
method computationally very time consuming. However, a revolving aver- 
aging procedure used in our previous work 4 seems to alleviate the need for 
using a large number of stochastic particles, N m , in any one given time step. 
In this averaging scheme, the solution provided to the CFD solver is based on 
an average of all the particles present over the last N av time steps instead of 
an average based solely on the number of particles present in any one single 
time step. This approach seemed to provide smooth Monte Carlo solutions 
to the CFD solver and, thereby, improving the convergence of the coupled 
CFD and Monte Carlo computations. The reason for improvement could be 
attributed to the average being based on N av N m particles instead of N m . 
Here, it is assumed that the solution contained in different time steps to be 
statistically independent of each other. 

V. Parallelization 

There are several issues associated with the parallelization of the PDF 
computations. The goal of the parallel implementation is to extract maxi- 
mum parallelism so as to minimize the execution time for a given application 
on a specified number of processors. 16 Several types of overhead costs are as- 
sociated with parallel implementation which include data dependency, com- 
munication, load imbalance, arithmetic, and memory overheads. Arithmetic 
overhead is referred to the extra arithmetic operations required by the paral- 
lel implementation and memory overhead refers to the extra memory needed. 
Excessive memory overhead reduces the size of a problem that can run on a 
given system and the other overheads result in performance degradation. 16 
Any given application usually consists of several different phases that must 
be performed in certain sequential order. The degree of parallelism and data 
dependencies associated with each of the subtasks can vary widely. 16 The 
goal is to achieve maximum efficiency with a reasonable programming effort. 

In our earlier work, we discussed the parallel implementation of a spray 
algorithm developed for the structured grid calculations on a Cray T3D. 6 
These computations were performed in conjunction with the application of 
the Monte Carlo PDF method to spray flames. The parallel algorithm made 
use of the shared memory constructs exclusive to Cray MPP (Massively Par- 
allel Processing) Fortran and the computations showed a reasonable degree 
of parallel performance when they were performed on a NASA LeRC Cray 
T3D with the number of processors ranging between 8 to 32. Later on, the 
extension of this method to unstructured grids and parallel computing in 
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PE 2 


PE 3 


[ 4 — PE 4 ^1 


Fig. 1 An illustration of the parallelization 
strategy employed in EUPDF. 


Fortran 77 with PVM calls was reported in Ref. 7. The Fortran 77 ver- 
sion offers greater computer platform independence. In this section, we only 
highlight some important aspects of parallelization. 

In the domain decomposition employed, the domain of computation is 
simply divided into n-parts of equal size and each part is solved by a different 
processor. Fig. 1 illustrates a simple example of the domain decomposition 
strategy adopted for the gas-phase computations where the total domain is 
simply divided equally amongst the available computer processing elements 
(PEs). In this case, we assumed the number of available PEs to be equal to 
four. 

At the beginning of the computation, all the information that is needed 
from the adjoining cells in computing the cell-face variables at the boundaries 
of the interface separating two neighboring processor-domains is obtained 
from the appropriate processors and stored in ghost cells. And interproces- 
sor communications are performed only at the beginning of integration step 
as the need arises. With the domain decomposition adopted, all the stages of 
a single PDF integration step involving the spatial transport, molecular mix- 
ing? s P r ay? and chemical kinetics lend themselves perfectly to parallel com- 
puting. Therefore, the Monte Carlo simulation is ideally suited for parallel 
computing and the run time could be considerably minimized by performing 
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START 


SPAWN MULTIPLE PEs FOR PARALLEL COMPUTING 

— i z 

READ PARAMETERS & GEOMETRIC DATA 


CALL SPRAY_PDF_READ PARAMETERS 


INITIALIZE r^^ A _ L L P _PL J ^?_ E ^J 
[INITIALIZE PDF VARIABLES J |READP_ 


Listed in Appendix VI 


RERUN 


READ PDF RESTART FILES 


INITIALIZE 


CALL SPRAY INT RERUN 


[RERUN 


INITIALIZE SPRAY VARIABLES I READ SPRAY RESTART FILES 


RERUN 


INITIALIZE 


RERUN 


INITIALIZE FLOW VARIABLES 



READ FLOW RESTART FILES 


UPDATE THERMODYNAMIC & TRANSPORT PROPERTIES 
r caUldclr to"advance SP^AYlQUATOllS 1 

L — — — — — — — — —^nzr — — — — — — — — — J 

I CALL PDF TO ADVANCE PDF EQUATIONS'! 

L , 1 

ADVANCE FLOW EQUATIONS I 



*CALL SPRAY PDF OUTPUT _ r^ Llsted in A PP endlx VI1 


OUTPUT 


STOP 


Fig. 2 The overall flow structure of the combined flow, LSPRAY, ad EUPDF solvers. 
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the computations on a massively parallel computer. 

VI. Details of the Coupling With the Flow &; Spray Solvers 

The PDF module is designed so that it could easily be coupled with 
any of the existing unstructured-grid flow and spray solvers. If geometric 
grid parameters - e.g. area vectors, grid connectivity, etc., were supplied 
separately, it could even be coupled with any of the existing structured-grid 
flow solvers. However, the present release of the code relies on the other 
modules of NCC for obtaining that information. 

The structure of the PDF solver is so designed that only a minimal 
amount of code modifications needs to be made within the flow and spray 
solvers for their coupling with EUPDF. The present version of the module 
relies entirely on the use of the Fortran common blocks for information ex- 
change between the various modules. Even this reliance should entail only 
few changes to be made within the PDF code for linkage with different solvers. 

The coupling issues could be understood better through the use of a 
flow chart shown in Fig. 2. The chart contains several blocks - some shown 
in solid lines and others in dashed lines. The ones in solid blocks represent 
the flow chart that is typical of most flow solvers. The ones in dashed blocks 
represent the coupling required for adding the spray and PDF solvers. The 
details on the spray blocks are not provided in this report as they could 
be found in the separate reference. 11 It should be borne in mind that the 
PDF solver could be run without the spray solver and vice-versa as they are 
independent. 

The flow chart for a typical flow solver begins by calling several routines 
- some for initiating the established PVM protocol for parallel computing 
and the others for spawning children of the same processes so that the com- 
putations could be performed simultaneously on various PEs participating 
in the parallel computing environment. It is followed by a routine to read 
various initial parameters. The geometric data could be either read directly 
or created by the inclusion of appropriate calling routines needed for grid 
generation. Then, the initial conditions for the flow variables need to be 
either specified or read from the restart files if it is a rerun. The thermody- 
namic and transport properties are then updated before advancing the flow 
equations over a series of time steps until the desired number of iterations 
are reached. Finally, the program is terminated after writing the output data 
on the separate restart and standard files. 

The coupling starts with the addition of a calling routine - 
spray _pdf _read_parameters - to read the spray and PDF control parame- 
ters followed by calls to the restart or initialization routines: pdf_int_rerun 
followed by spray _int_rerun. Then, calls to dclr and pdf were made in 
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order to advance the spray and PDF equations in a sequential order be- 
fore advancing the flow equations. It should be borne in mind that if the 
PDF solver is invoked, the thermodynamic and transport properties would 
be evaluated by the routines contained within the PDF solver instead of the 
ones contained in the flow solver. Finally, spray _pdf_output and outpdf2 
routines are included for outputting the PDF and spray data on appropriate 
restart and standard files. 

Appendix I contains a listing of all of the Fortran subroutines and func- 
tions used in EUPDF. 

A brief description of each one of these routines is provided in Appendix 

II. 

Appendix III contains the listing of a subroutine which is used for read- 
ing some of the control and other associated parameters involving LSPRAY 
and EUPDF solvers. 

Appendix IV contains a list of some Fortran variables used by EUPDF. 

Appendix V contains a list of the geometric variables used by EUPDF 
which are currently supplied by the flow code of NCC. 

Appendix VI contains a subroutine for initialization and restart of the 
PDF computations. For the case of initialization, the PDF variables are 
initialized based on the solution provided by the flow solver. For the case of 
restart, the the PDF variables are read from a previous solution. 

Appendix VII contains the listing of subroutines spray _pdf_output 
and outpdf2 which are for writing output data from LSPRAY and EUPDF 
codes on separate standard and restart files. 

Appendix VIII contains an example of the partial listings of code initi- 
ation for coupling LSPRAY and EUPDF with a gas flow solver. 

Appendix IX contains a listing of all subroutines and functions used 
in the evaluation of the thermodynamic and transport Properties which are 
used not only in EUPDF but also in the flow solver as well as in LSPRAY. 
These subroutines include rctinp, props_ther, props_tran, 
get_kt_and_cp_loc_mcs, find_inverse_molecular_weight_mcs, and 
find_h_mcs. 

The last appendix provides an example of the summary of the CPU times 
taken by CORSAIR and EUPDF for the case of a confined swirl-stabilized 
spray flame when the computations were performed on a LACE cluster at 
NASA LeRC. 

VII. Concluding Remarks 

A Monte Carlo technique has been outlined for the computation of spray 
flames on unstructured grids with parallel computing. The method outlines 
several techniques designed to overcome some of the high computer time and 
storage limitations associated with the Monte Carlo simulation of practical 
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combustor flows. The viability of the present method for its application to 
the modeling of practical combustion devices is demonstrated through the 
ease with which grids could be generated for complex combustor geometries 
by the use of automated unstructured mesh generators and the ability to run 
the computations on massively parallel computers. 
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Appendix I 


Table 1. A List of EUPDF Fortran Subroutines and Functions. 

Number 

Name of the routine 

Code released with 

Code page/ 



NCC/Description page 

User’s manual 

i. 

chem2 

25 


2. 

coeffivQ 

25 


3. 

convec 

25 


4. 

find_h_mcs() 

25 

61 

5. 

find Jnverse -molecular .weight _mcs() 

26 

60 

6. 

get_kt.and_cp_loc_mcs() 

26 

58 

7. 

molmix 

26 


8. 

outpdflQ 

26 


9. 

outpdf‘2() 

26 

46 

10. 

pdf 

26 


11. 

pdf_int_rerun 

27 

39 

12. 

prntQ 

27 


13. 

props.ther 

27 

55 

14. 

props.tran 

27 

57 

15. 

rctinp 

27 

53 

16. 

rep lac 

27 


17. 

spray 

27 


18. 

spray _pdf .output 

27 

43 

19. 

spray _pdf_read_parameters 

27 

29 
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Appendix II 


Description of EUPDF Fortran Subroutines and Functions 

1. subroutine chem2: 

PURPOSE: This subroutine integrates the chemical kinetics step 
of the PDF solution which is especially tailored 
to solve a single step global mechanism of the Westbrook 
and Dryer for n-heptane oxidation. 

2. subroutine coeffiv(axyzp,d): 

PURPOSE: Computes convection/diffusion coefficients based on 

a cell-centered finite-volume derivative for diffusion 
and a upwind or hybrid differencing scheme for the 
convection . 

3. subroutine convec: 

PURPOSE: This is the controlling routine for the integration 
of the convection/diffusion step. 

(1) Computes the convection/diffusion coefficients by 
calling the subroutine coeffiv. 

(2) Converts the decimal numbers into integer numbers 
which represent the number of particles to be 
transferred across the neighboring cells. 

(3) Loads the random numbers into appropriate arrays 
which are used in determining the particles that 
need to be transferred across neighboring cells. 

(4) Moves particles across different neighboring cells 
based on the computed convection/diffusion 
coefficients by calling the subroutine replac. 


4. function find_h_mcs(element, centroid): 

PURPOSE: Computes specific enthalpy at a nodal point of the 
computational grid. 
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5. function find_inverse_molecular_ weight_mcs(element, centroid): 

PURPOSE: Computes inverse of the molecular weight of a gaseous 
mixture at a nodal point of the computational grid. 

6. subroutine get_kt_and_cpJoc_mcs(ijk, centroid): 

PURPOSE: Computes specific heat and thermal conductivity of a 
gaseous mixture at a nodal point of the computational 
grid. 

7. subroutine molmix: 

PURPOSE: This is a routine for computing transport in scalar space 
due to molecular mixing. It uses the relaxation to the 
ensemble-mean molecular mixing model. 

8. subroutine outpdfl(ncyc): 

PURPOSE: Calculates the residuals for the PDF solution. 

9. subroutine outpdf2(ncyc): 

PURPOSE: Outputs the PDF solution to restart files. 

10. subroutine pdf: 

PURPOSE: This is the main controlling routine for the 
Monte Carlo PDF solver. It updates the pdf 
solution through a series of calls to different 
subroutines and returns control back to the 
calling routine. It also performs the following 
functions : 

(1) controls interprocessor communications. 

(2) calculates the average value of scalars over particles. 

(3) transfers the time-averaged PDF solution into 
corresponding arrays of the conventional CFD 
solver . 

(4) computes thermodynamic and transport properties. 

(5) computes the residual terms. 
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11. subroutine pdfJnt_rerun: 

PURPOSE: This routine is used either to initialize the Monte 

Carlo PDF computations or to restart the PDF computations 
from the output of a previous PDF calculation. 

12. subroutine prnt(aass,phig,nbc,nnode,np,ipid, nstart,nend,iul): 

PURPOSE: I/O print out. 

13. subroutine props_ther: 

PURPOSE: This routine computes enthalpy, specific heat, and 
density . 

14. subroutine props.tran: 

PURPOSE: This routine computes transport properties. 

15. subroutine rctinp: 

PURPOSE: Initializes parameters used in the thermodynamic 
and transport property evaluation as well as in 
the chemical kinetics scheme. 

16. subroutine replac: 

PURPOSE: Moves particles across different neighboring cells 
during the convection/diffusion step. 

17. subroutine spray: 

PURPOSE: Updates PDF solution associated with the contribution 
arising from the liquid-phase source terms. 

18. subroutine spray _pdf_output: 

PURPOSE: This routine writes output data from EUPDF & LSPRAY 
computations to restart and standard-output 
files . 

19. subroutine spray _pdfLread_parameters: 

PURPOSE: This routine reads controlling parameters associated 

with the EUPDF and LSPRAY solvers. Based on the controlling 
parameters read, it might invoke an initialization routine 
of the EUPDF solver which is needed in the thermodynamic Sc 
transport properties evaluation. 
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Appendix III 

A Subroutine Listing For the Read Parameters of LSPRAY and 

EUPDF 


c 

subroutine spray _pdf _read_parameters 
c 

include * dcf slog . i J 
include ’ dcf slog_rw . i * 
c 

c 

c 

c 

c PURPOSE: This routine reads controlling parameters associated 
c with the EUPDF and LSPRAY solvers. Based on the controlling 

c parameters read, it might invoke an initialization routine 

c of the EUPDF solver which is needed in the thermodynamic Sc 

c transport properties evaluation, 

c 

c FORM OF CALL: call spray_pdf _read_parameters 

c 

c 

c ADDITIONAL I/O: 
c 

c INPUT: spray_pdf _parameter_input 

c 

c OUTPUT : None 

c 

c 

c 

c lspray controls turning on or off spray computations, 
c lspray = .TRUE. - turns on spray computations, 
c = .FALSE. - otherwise, 

c 

c ldread controls reading or not from restart files for 

c spray computations. 

c ldread = .TRUE. - restarts from previous runs, 
c = .FALSE. - starts from initial conditions, 

c 

c ispray_mod= This variable controls calls to the spray 
c solver. The spray solver is called once at 

c every ispray_mod times of CFD iterations. 
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c 

c ipread = Assigned unit number for the file: liquid. input . 

c idread = Assigned unit number for the file: liquid.results . 

c idread2= Assigned unit number for the file: liquid_results.ini 

c idwrit = Assigned unit number for the file: liquid_results_new 

c idwrit2= Assigned unit number for the file: liquid.results.ini 

c 

c ipdf controls turning on or off Monte Carlo PDF computations, 
c ipdf = 0 turns off Monte Carlo PDF computations, 
c =1 otherwise, 

c 

c ns serves two functions depending on whether ns has a 
c zero or non-zero value. 

c ns = 0 starts the Monte Carlo PDF computations from 
c initial conditions. 

c ns = a non-zero number restarts the computations from 

c a previous run. a non-zero number represents the 

c last iteration number of a previous run which is 

c used in the time-averaging scheme utilized in 

c the PDF computations, 

c 

c ipdf.mod = This variable controls calls to the PDF 
c solver. The PDF solver is called once at 

c every ipdf.mod times of CFD iterations, 

c 

c ipdf_num = In a given cycle, the pdf solver is advanced over 
c a number of time steps given by ipdf.num. 

c 

c ireal = Assigned unit number for the file: pdf .results, 

c irea2 - Assigned unit number for the file: pdf .results.ave . 

c iwril = Assigned unit number for the file: pdf .results, 

c iwri2 = Assigned unit number for the file: pdf .result s.ave . 

c 

c 

c 

open(unit=85 ,f ile= J spray.pdf .parameter. input ; ) 
read(85 , *) 

read(85 , *) 1 spray , Idread, ispray.mod 
read(85,*) 

read (85 , *) ipread , idread , idwrit , idread2 , idwrit2 
read(85 ,*) 

read (85 ,*) ipdf , ns , ipdf .mod, ipdf _num 
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c 


read(85 , *) 

read (85 ,*) ireal , irea2 , iwril , iwri2 
close(unit=85) 


c 

c 

c 

c Routine rctinp of the PDF solver provides initialization 
c parameters used in the themodynamic and transport property 
c evaluation as well as in the chemical kinetics scheme, 
c 

if (ipdf . eq . 1) then 
call rctinp 
endif 
c 

c 

c 

RETURN 

END 

c 


31 




Appendix IV 


A Listing of Some Fortran Variables Used in EUPDF 

c 

c 

c 

c arh = A = constants in the equation used for representing the 
c chemical reaction rate, 

c 

c avl() It contains an average based on nparti*nav particles, 

c It is an average employed in the revolving time-averaging 

c technique, 

c 

c avmanyO It contains nav number of individual averages which 
c are computed over nparti. 

c 

c centroid = logical variable. If yes, it repesents the values at 
c the cell center. If not, it represents values at the 

c cell faces . 

c 

c cp0,..cp4 = coefficients of the polynomial used in determining the 

c variable specific heat, 

c 

c dhform = coefficients of the polynomial used in determining the 

c specific enthalpy, 

c 

c dwm = coefficients of the polynomial used in determining the 

c inverse of the molecular weight of a mixture, 

c 

c eno2, epox, 
c epfu,emno2, 

c empox,empfu = constants used in the chemical 
c kinetics solver, 

c 

c erh = E_a = constants in the equation used for representing the 
c chemical reaction rate, 

c 

c f 1 ( ) = fuel mass fraction at the center of the grid cell, 

c flbarO = fuel mass fraction at the cell face, 
c 

c f2() = oxygen mass fraction at the center of the grid cell. 
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c f2bar() * oxygen mass fraction at the cell face. 

c 

c 

c fluidO = logical variable. If yes, it repesents a fluid 
c element. If not, it represents a solid element, 

c 

c h() = enthalpy at the center of the grid cell, 

c hbar() = enthalpy at the cell face, 
c 

c iconpc = store differences in array index convention used by the 
c flow and pdf codes, 

c 

c nav = N_{av} = number of time steps employed in the time-averaging technique. 

c 

c 

c nparti = N_m = number of particles per cell 

c 

c 

c ns = current iteration of the PDF solver, 

c 

c nscala = sigma = number of scalars 
c 

c parti contains pdf solution from the previous time step 
c 

c part2 contains pdf solution for the current time step 
c 

c pz = pressure 
c 

c ru = universal gas constant 
c 

c schmdi= Sc = Schmdit number (=1/0.7). 
c 


c t() = temperature at the center of the grid cell, 

c tbar() = temperature at the cell face, 
c 

c wl,w2,w3 = constants resulting from some algebraic manipulation 
c which reduces the need for solving speciess equations 

c from 5 to 3. 

c 

c wmole = W = molecular weight 


c xnup, xnupp = nu = stoichiometric coefficients of reactants and 
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Appendix V 

A Listing of Geometric Variables Used in LSPRAY and EUPDF 


c 

c 

c 

c The EUPDF and LSPRAY modules expects the following inputs on 
c the grid related information: 
c 

c nodes = total number of the computational elements, 

c nedge = total number of cell faces in the computational domain, 

c nfaces(i) = total number of faces of the element, 1 . 

c 

c edge(i,l) and edge(i,2) represent the adjacent elements 
c of the face, i, if the face, i, if the face happens 

c to be an interface between two elements. Otherwise 

c edge(i,l) represents the correponding boundary 

c condition identifier, 

c 

c f ace_to_edge(i , j) represents the face ID of the element, i, 
c and the face , j . 

c 

c cl(i,j) provides connect ivity map . cl(i,j) = adjacent element 
c ID of the element, i, and the face, j, otherwise cl(i,j) 

c = boundary condition identifier on any boudary . 

c 

c vol(i) = volume of the element, i. 
c 

c areax(i) , areay(i), and areaz(i) are the cartesian components 
c of the outward pointing area vector of the face, i. 

c 

c xl(i), yl(i), zl(i) are the cartesian components of the node 
c one of the element, i. Similarly, x2(i), y2(i), z2(i) 

c are for node 2 and so on. 

c 

c triangle(i) is .true, if i is a triangular element. Similarly, 
c quadrilateral (i) , tetrahedron(i) , and wedge(i) are logical 

c varibles representing other type of elements, 

c 

c axisymmetric is set to .true, for axisymmetric computations 

otherwise it is .false. The axisymmetric computations are 
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performed by generating 3D elements from a 2D mesh with 
an arc centered around the z coordinate, z=0.0. The angle 
of the arc is defined by the variables, ARC, in radians 
and THETAO, in degrees. 



Appendix VI 

A Subroutine Listing For EUPDF Code Initialization and Restart 


c 

subroutine pdf _int_rerun 
c 

include ' cf sparms . i } 
include ' cfsdt.i' 
include > cf spert . i * 
include * cf sconv . i ; 
include } cf stime . i * 
include * cf smimd . i * 

include * cf sarea . i ; 
include ; cf snodes . i J 
include } dcf slog * 1 J 
include } dcf slog_rw . i } 
include ' cf svars . i * 
include ; cf sprop . i * 
include ; cfsh.i J 
c 

c Include common blocks associated with PDF computations, 
c 

include 'p3dpar . i ; 
include ; p3dcom. i ; 
include ; p3dave . i > 
include ; p3dpro . i ' 
c 

c 

c 

c 

c PURPOSE: This routine is used either to initialize the Monte 
c Carlo PDF computations or to restart the PDF computations 

c from the output of a previous PDF calculation, 

c 

c FORM OF CALL: call pdf _int_rerun 
c 

c 

c 

c 

c 

c Initialize the computations if it is not a rerun. 
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c 

c Initialize the PDF particle attributes as well as the 
c averages to those provided by the flow solver. 

c 

if (ipdf .eq.l.and.ns.eq.O) then 

do 39 ijk = 1, nodes 

avyl=f 2(ijk) 

avy2=f 1 (ijk) 

avye=h(i jk) 

avyt=t (ijk) 

avmany(ijk,pox,2) = avyl 
avmany (ijk, pfu, 2) = avy2 
avmany (ijk, pen, 2) = avye 
avmany (ijk,pte ,2) = avyt 
avmany(ijk,pox, 1) = avyl 
avmany ( i jk, pf u, 1) = avy2 
avmany (ijk, pen, 1) = avye 
avmany (i jk, pte, 1) = avyt 
39 continue 

do 391 ijk = 1, nodes 
do 391 ip = l,nparti 
partl(ijk,ip,pox) = avmany(ijk,pox,2) 
parti (ijk, ip, pfu) = avmany(ijk,pfu,2) 
parti (ijk, ip, pen) = avmany(ijk,pen,2) 
parti (ijk, ip, pte) = avmany (ijk, pte, 2) 
part2(ijk, ip, pox) = avmany (ijk, pox, 2) 
part2(ijk,ip,pfu) = avmany (ijk, pfu, 2) 
part2(ijk, ip, pen) = avmany (ijk, pen, 2) 
part2 (i jk, ip, pte) = avmany (ijk, pte, 2) 

391 continue 

do 321 ijk=l, nodes 
do 321 n=l,nav 

avl(ijk,pox,n) = avmany (ijk, pox, 2) 
avl (ijk, pfu, n) = avmany(ijk,pfu,2) 
avl(ijk,pen,n) = avmany (ijk, pen, 2) 
avl (ijk, pte ,n) = avmany (ijk, pte, 2) 

321 continue 
ns=nav 
c 

c Compute thermodynamic properties, 
c 

call props_ther 
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c 

c Compute transport properties, 
c 

call props_tran 
endif 
c 

c 

c 

c 

c Read restart files. Also, initialize some other parameters 
c of interest, 
c 

if (ipdf . eq. 1 . and .ns .gt .nav) then 
itrecl=nparti*nscala*4+4 
open(unit=ireal ,f ile= J pdf .results y , 

>access= > direct * , recl=itrecl ,f orm= Unformatted > ) 
itrecl=nscala*nav*4+nscala*2*4+4 
open(unit=irea2,f ile= 'pdf .result s_ave 7 , 

>access= * direct * , recl=itrecl ,f orm= } unformatted ' ) 
do ijk=l, nodes 
irecord=ijk+nodes* (ipid-1) 
read(ireal , rec=irecord) 

> ( (part 2 (ijk , il,i2) , i 1 = 1 ,npart i) , i2=l ,nscala) 
enddo 

do ijk=l, nodes 
irecord=ijk+nodes* (ipid-1) 
read(irea2 ,rec=irecord) 

>((avl(ijk,il,i2) , i 1=1 , ns cal a) , i2=l ,nav) , 

>( (avmany (ijk, i3, i4) , i3=l , ns cal a) , i4=l , 2) 
enddo 
c 

do 4000 ijk=l, nodes 
do 4000 i=l,nparti 
do 4000 k3=l,nscala 
4000 part 1 (ijk, i ,k3)=part2 (ijk, i ,k3) 
c 

c Compute thermodynamic properties, 
c 

call props_ther 
c 

c Compute transport properties, 
c 
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c 

c 

c 

c 

c 


c 


call props_tran 
endif 


return 

end 
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Appendix VII 


A Subroutine Listing For LSPRAY and EUPDF Data Output 


subroutine spray_pdf .output 
c 

include 9 cf sparms . i 9 
include ; cfsdt.i J 
include 9 cf spert . i ; 
include 9 cf sconv . i 9 
include } cf st ime . i * 
include } cf smimd . i ' 

include 9 cf sarea . i J 
include * cf snodes . i 9 
include 9 cf svars . i ' 
include ; cf sprop . i ; 
include ' cfsh.i ; 
c 

c Include common blocks associated with spray and PDF computations, 
c 

include 9 dcf slog . i 9 
include ' dcf slog.rw . i * 
c 

c Include common blocks associated with PDF computations, 
c 

include * p3dpar . i } 
include 9 p3dcom . i ' 
include ; p3dave . i ; 
include 9 p3dpro . i ’ 
c 

c Include common blocks associated with spray computations, 
c 

include ' d3dpar . i ; 
include 9 d3dcom . i } 
include ’ d3dinj . i ; 
include * d3dprl . i ; 
c 

c 

c 

c 
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c PURPOSE: This routine writes output data from PDF & spray 
c computations on a separate restart and standard 

c output files, 

c 

c FORM OF CALL: call spray_pdf _output 

c 

c 

c ADDITIONAL I/O: 
c 

c INPUTS : None . 

c 

c OUTPUTS : 

c 

c liquid_results_new 

c liquid_results_ini 

c spray_pdf .parameter. input 

c 



c 

c 

c 

c 

c Write spray restart files . 
c 

if(lspray) then 

open(unit=idwrit ,f ile=' liquid.results.new’ , 

> access= ' direct ’ ,recl=136 , 

> f orm= 'unf ormatted' ) 
if (ipid.eq. 1) then 

open(unit=idwrit2,f ile =, liquid_results_ini ’ ) 
write(idwrit2,*)nr_total 
call f lush(idwrit2) 
write(idwrit2,*)dtil,dtml,tl ,tll,tml 
call f lush(idwrit2) 
write(idwrit2 , *) iseed 
call f lush(idwrit2) 
close (unit=idwrit2) 
endif 
c 

INTS_DATA_2=314 
do n=l,np 
no_to_ip(n)=nr 
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enddo 

do n=l,np 

if (ipid.ne.n) then 

irc= send_data_i (iul (n) ,no_to_ip(n) , 1, INTS_DATA_2) 

endif 

enddo 

do n=l,np 

if (ipid.ne.n) then 

irc= recv_data_i (iul (n) ,no_f r_ip(n) , 1, INTS_DATA_2) 

endif 

enddo 

no_f r_ip(ipid)=no_to_ip(ipid) 

irecordd=0 

do n=l , ipid-1 

irecordd=irecordd+no_f r_ip (n) 
enddo 
c 

do ip=l,nr 
irecord=irecordd+ip 
isent=isen(ip) +(isep(ip) -1) *nodes 
write(idwrit ,rec=irecord) ndrr(ip) ,ins(ip) , 

1 isent ,xki(ip) , yki(ip) ,zki(ip) ,uki(ip) , 

2 vki(ip) ,wki(ip) ,tki(ip) ,rki(ip) ,ski(ip) ,sklim(ip) , 

3 (vh(ip , j ) , j = 1 , nde+4) 
call f lush(idwrit) 

if(ip.ge.l) then 

write (1 ,*) irecord,ndrr (ip) , ins (ip) , 

1 isent ,xki(ip) ,yki(ip) ,zki(ip) ,uki(ip) , 

2 vki(ip) ,wki(ip) ,tki(ip) , rki(ip) , ski (ip) , sklim(ip) , 

3 (vh(ip,j) ,j = l, nde+4) ,nr ,nr_total , irecord 
endif 

enddo 

close (unit=idwrit) 
endif 
c 

c 

c 

c 

c 

c 

c Update file: spray_pdf .parameter^ input . 
c Also, write PDF restart files. 
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c 

c 

if (ipdf .eq. 1) then 
if (ipid.eq. 1) then 

open(unit=85 ,f ile= 7 spray_pdf _parameter_input 7 ) 

write(85,*) 7 lspray ldread ispray_mod 7 

write (85 , *) lspray , ldread , ispray_mod 

write (85,*) 7 ipread idread idwrit idread2 idwrit 2 7 

write (85 , *) ipread , idread , idwrit , idread2 , idwrit 2 

write (85 ,*) 7 ipdf ns ipdf_mod ipdf_num 7 

write (85 , *) ipdf , ns , ipdf _mod , ipdf _num 

write (85 , *) 7 ireal irea2 iwril iwri2 7 

write (85 , *) ireal , irea2 , iwril , iwri2 

close(unit=85) 

endif 

call outpdf2(ns) 
endif 
c 



c 



c 

c 

c Write output of spray computations either to unit 
c one or to the screen, 
c 

if (lspray) call prnspr 
c 

return 

END 

c 

subroutine outpdf 2 (ncyc) 
include 7 p3dpar . i 7 

include 7 cf sparms . i 7 
include 7 cf snodes . i 7 


include 7 p3dcom . i 7 
include 7 p3dave . i 7 


include 7 parallel . i 7 
include 7 cf smimd . i 7 
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include } dcf slog_rw . i ; 


c 

c 

c 

c PURPOSE: Outputs the PDF solution to restart files, 
c 

c FORM OF CALL: call outpdf2() 
c 

c ARGUMENTS/PARAMETERS: ncyc 
c 

c ADDITIONAL I/O: pdf. results, pdf .result s.ave 
c 

c OUTPUT: outputs PDF solution to pdf .results & 
c outputs PDF averaging-procedure solution 

c to pdf .results.ave . 

c 

c 

itrecl=nparti*nscala*4+4 
open(unit=iwril , f ile= ; pdf .results ' , 
>access='direct J ,recl=itrecl ,form= 'unformatted J ) 
c 

itrecl=nscala*nav*4+nscala*2*4+4 
open(unit=iwn2,f ile= 'pdf .results.ave' , 

>access= 'direct ' ,recl=itrecl ,f orm= 'unformatted' ) 
c 

do ijk=l, nodes 
irecord=ijk+nodes* (ipid-l) 
write(iwril ,rec=irecord) 

> ( (part 2 (ijk, il , i2) , i 1= 1 ,nparti) , i2=l ,nscala) 
call flush(iwril) 
enddo 
c 

do ijk=l, nodes 
irecord=ijk+nodes*( ipid-l) 
write(iwri2 ,rec=irecord) 

>((avl(ijk,il,i2) , il = l , ns cal a) , i2=l ,nav) , 

> ( (avmany (i jk, i3 , i4) , i3=l ,nscala) , i4=l ,2) 
call flush(iwri2) 
enddo 
c 
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close (unit=iwril) 
close (unit=iwri2) 


c 

return 

end 

c 

c 



Appendix VIII 


An Example of the Partial Listings of Code Initiation For 
Coupling LSPRAY and EUPDF With a Gas Flow Solver 

1. The following segment shows how include calls to 
spray _int_rerun &; pdf_int_rerun. 


c 

c 

c 

c 

c Include common blocks associated with spray and PDF computations, 
c 

include , dcf slog . i ' 
include 9 dcf slog.rw . i } 
c 
c 

c Initialize Monte Carlo PDF computations, 
c 

if (ipdf . eq . 1) then 
call pdf _int_rerun 
endif 
c 

c Initialize spray computations, 
c 

IF(lspray) then 
call spray_int_rerun 
endif 
c 

c 

II. The following segment shows how to include calls to DCLR &; 
PDF. 


c 

c Include common blocks associated with spray and PDF computations, 
c 

include 9 dcf slog . i 9 
c 

double precision tbiggas, tendgas, totaltgas 
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c 


c 

c 

c 

c Call dclr in order to advance the spray computations 
c over a time step of dtgl . 
c 

if (lspray.and.mod(iteration,ispray_mod) .eq.O) then 
call dclr 
endif 
c 



c 



c 

c 

c 

c Call pdf in order to advance the PDF computations 
c over the next time step, 
c 


c 

c 

c 


if (ipdf . eq . 1 . and . mod (iter at ion, ipdf _mod) .eq.O) then 
do i=l,ipdf_num 
call pdf 
enddo 
endif 


III. The following segment shows how to include the interphase 
contributions to the gas phase governing equations. 




c 

c 

c Include common blocks associated with spray and PDF computations, 
c 

include ' dcf slog . i } 
c 

c Include common blocks associated with spray the solver, 
c 
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include ; d3dqat . i * 
c 

c 

c 

c 

c 

c 

c Include liquid-phase contributions to mass, momentum, species, 
c and energy equations, 
c 

if(lspray) then 
do i=l, nodes 

sourcem(i) =sourcem(i) +qmsx(i) 
sourceu(i) =sourceu(i)+qmsx(i) 
sourcev(i)=sourcev(i)+qmsy (i) 
sourcew(i)=sourcew(i)+qmsz(i) 
sourcef ( i ) =sourcef ( i ) +qms ( i ) 
sourceh ( i ) =sourceh (i ) +qmsh (i ) 
enddo 
endif 
c 

c 
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Appendix IX 


A Listing of Subroutines Used For Thermodynamic and 
Transport Property Evaluation 


c 

c 

subroutine rctmp 

include J p3dpar . i * 

include } cf sparms . i J 
include * cf snodes . i J 

include ’ p3dcom . i J 
include * p3dpro . i * 
c include ; p3dtim.i ; 

dimension cptO(nspt ,3) ,cptl (nspt ,3) , cpt2 (nspt , 3) , 

> cpt3(nspt,3) ,cpt4(nspt ,3) 

c 

c 

c 

c PURPOSE: Initializes parameters used in the themodynamic 
c and transport property evaluation as well as in 

c the chemical kinetics scheme. 

c 

c 

c ru = universal gas constant 
c pz = pressure 
c 

ru =8314.3 

pz =1.01325e+5 

cppdf =1297.1012 
tref pdf =298 . 0 
c 

c default species: 


c 



c 

o2 

= 1 

c 

c7hl6 

= 2 

c 

h2o 

= 3 

c 

co2 

= 4 

c 

n2 

= 5 
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c 

c iconpc = store differences in array index convention used by the 
c flow and pdf codes . 

c 

iconpc(l )=2 
iconpc(2)=l 
iconpc(3)=3 
c 

c read datafile on properties for n-heptane. 
c 

c wl,w2,w3 = constants resulting from some algebraic manipulation 
c which reduces the need for solving speciess equations 

c from 5 to 3. 

c 

c xnup, xnupp = stoichiometric coefficients of reactants and 
c products, 

c 

c arh,erh = constants in the equation used for representing the 
c chemical reaction rate, 

c 

c wmole = molecular weight 

c 

c cp0,..cp4 = coefficients of the polynomial used in determining the 

c variable specific heat, 

c 

c dwm = coefficients of the polynomial used in determining the 

c inverse of the molecular weight of a mixture, 

c 

c dhform = coefficients of the polynomial used in determining the 

c specific enthalpy, 

c 

c eno2, epox, 
c epfu,emno2, 

c empox,empfu = constants used in the chemical 
c kinetics solver, 

c 

open (unit =35 ,file =, nheptane_property_file ' ) 

read (35 , *)wl , w2, w3 

read (35,*) (xnup(l ,i) ,i=l,4) 

read (35,*) (xnupp (l,i) ,i=l,4) 

read (35 , *) arh(l) ,erh(l) 
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read(35,*) (wmole(i) ,i=l,5) 
c 

read(35 , *) C(cpO(i, j) , i=l ,5) , j=l ,3) , 

1 ((cpl (i , j ) , i=l ,5) , j=l ,3) , 

1 ((cp2(i,j) ,i=l,5) ,j=l,3) , 

1 ( (cp3(i ,j),i=l,5),j=l,3), 

1 ( (cp4(i ,j),i=l,5),j=l,3) 

c 

read (35 , *) (dwm(i) , i=l ,3) 
c 

read(35,*) (dhform(i) , i=l ,3) 
c 

read (35 , *) eno2 , epox , epf u , emno2 , erapox , empf u 
close(unit=35) 
c 
c 
c 

return 

end 

subroutine props_ther 
include ' cf sparms . i ’ 
include ; cf snodes . i J 
include * cf svars . i } 
include J cfsh.i ; 
include ; cf sprop . i * 
include J cf spref . i ; 
include ; cf sgas . i ' 
include > cf scpmix . i ' 
include ; con j ugate . i * 
c 

include y p3dpar . i * 
include ; p3dpro . i } 

c 

c 

c PURPOSE: This routine computes enthalpy, specific heat, and 
c density. 

c 

c 

do 100 ijk=l, nodes 
if (f luid(ijk) ) then 
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tmp= t(ijk) 

icspec=iconpc(3) 

dO=dhf orm(icspec) 

icspec=iconpc(l) 

dO=dO+ f 1 (ijk)*dhform(icspec) 

icspec=iconpc(2) 

dO=dO+ f 2 (i jk) *dhf orm(icspec) 

do 1010 iter=l,l 

j c=cvmgp (2,1, tmp-600 . 0) 

jc=cvmgp(3, jc,tmp-1000 .0) 

icspec=iconpc(3) 

dl= cpO(icspec, jc) 

d2= cpl(icspec, jc) 

d3= cp2(icspec, jc) 

d4= cp3(icspec, jc) 

d5= cp4(icspec, jc) 

icspec=iconpc ( 1) 

dl=dl+ f l(ijk)*cpO(icspec,jc) 

d2=d2+ f l(ijk)*cpl(icspec, jc) 

d3=d3+ f l(ijk)*cp2(icspec, jc) 

d4=d4+ f l(ijk)*cp3(icspec, jc) 

d5=d5+ f 1 (ijk) *cp4(icspec , jc) 

icspec=iconpc(2) 

dl=dl+ f2(ijk)*cp0(icspec, jc) 

d2=d2+ f2(ijk)*cpl(icspec, jc) 

d3=d3+ f2(ijk)*cp2(icspec, jc) 

d4=d4+ f2(ijk)*cp3(icspec, jc) 

d5=d5+ f 2 (ijk) *cp4(icspec , jc) 

h(ijk)=((((0. 20*d5*tmp+0 . 25*d4) *tmp+d3/3 . 0) *tmp+ 

> 0 . 50*d2) *tmp+dl)*tmp+dO 
devm= ( ( (d5*tmp+d4) *tmp+ 

> d3)*tmp+d2)*tmp+dl 
1010 continue 

icspec=iconpc (3) 

rraw = dwm(icspec) 

icspec=iconpc(l) 

rraw = rmw + f 1 (ijk)*dwm(icspec) 

icspec=iconpc(2) 

rmw = rmw + f 2(ijk)*dwm(icspec) 

rho ( i j k) = (p (i j k) +p_ref erence) / (Rgas*trap*rmw) 

cp(ijk) = devm 

endif 
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100 continue 
return 
end 


subroutine props.tran 


c 


include 

include 

include 

include 

include 

include 

include 

include 

include 

include 

include 


* cf sparms . i * 

* cf snodes . i * 

, cf svars . i ; 

; cf smimd . i J 
; cf sprop . i ’ 

} cf sother . i J 
; cf sturb . i ' 

} cf svislam. i ' 
} cf sdt . 1 ' 

; cf scpmix . i J 
J conjugate . i ' 


c 

c 

c PURPOSE: This routine computes transport properties. 

c 

c 

c calculate the effective viscosity for turbulent flow 

if (turbulent) then 
lvis=0 


do 100 ijk=l, nodes 
if (fluid (ijk) ) then 

mu(ijk)=rho(ijk) *k(i jk) **2*cmu/ (e(ijk)+l . e-20) 

& +temp_vislam 

c enforce upper and lower bounds on turbulent viscosity 

if (mu (ijk) . gt . vismax) then 
lvis=l 

mu(ijk)=vismax 

endif 

mu(ijk)=amaxl (mu (ijk) , temp.vislam) 
kt (i jk)=mu(i jk)*cp(ijk) /prandl 
endif 

100 continue 

if (mod(itime, 50) . eq . 0 . and . ip id . eq . 1) then 
if (lvis . eq . 0) print 746 
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746 format (til , 'No limit enforced on turbulent viscosity . ; , /) 
if (lvis.eq. 1) print 747 

747 format (til , 'Limit enforced on turbulent viscosity . ; , /) 
endif 

else 

do 101 ijk=l, nodes 

if (fluid(ijk) ) then 

mu(ijk)=temp_vislam 

kt ( i j k) =mu ( i j k) *cp ( i j k) /prandl 

endif 

101 continue 
endif 
return 
end 

subroutine get^kt.and^cp.loc.mcsCijk, centroid) 
include } cf sparms . i ; 
include } cf snodes . i ; 
include * cf svars . i ; 
include } cf sprop . i J 
include ; cf sarea . i J 
include ; cf sf lux . i ; 
include * cf ssorc_enth . i J 
include ' cf schar . i ; 
include ' cf sadj . i > 
include , cfsdt.i > 
include ; cf sorder . i > 
include } cf scomp . i ; 
include ; cf scpmix . i } 
include ; conjugate . i * 
c 

include ' p3dpar . i } 
include 'p3dpro . i ; 

logical centroid 

c 

c 

c PURPOSE: Computes specific heat and thermal conductivity of a 
c gaseous mixture at a nodal point of the computational 

c grid. 

c 

c 
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if (centroid) then 
trap= t(ijk) 

jc=cvmgp(2 , 1 , tmp-600 . 0) 
j c=cvmgp (3 , j c , tmp-1000 . 0) 
icspec=iconpc (3) 
dl= cpO(icspec, jc) 
d2= cpl (icspec, jc) 
d3= cp2(icspec, jc) 
d4= cp3(icspec, jc) 
d5= cp4(icspec, jc) 
icspec=iconpc(l) 
dl=dl+ f 1 (ijk)*cpO(icspec, jc) 
d2=d2+ f l(ijk)*cpl(icspec, jc) 
d3=d3+ f 1 (ijk)*cp2(icspec, jc) 
d4=d4+ f 1 (ijk)*cp3(icspec, jc) 
d5=d5+ f 1 (i jk) *cp4(icspec , jc) 
icspec=iconpc(2) 
dl=dl+ f2(ijk)*cp0(icspec, jc) 
d2=d2+ f2(ijk)*cpl(icspec, jc) 
d3=d3+ f2(ijk)*cp2(icspec, jc) 
d4=d4+ f2(ijk)*cp3(icspec, jc) 
d5=d5+ f2(ijk)*cp4(icspec, jc) 
cp ( i j k) = ( ( (d5*tmp+d4) *tmp+d3) *tmp+ 
> d2)*tmp+dl 

kt (ijk)=mu(ijk) *cp(ijk) /prandl 
else 

tmp= tbar(ijk) 

j c=cvmgp (2 , 1 , tmp-600 . 0) 

jc=cvmgp(3 , jc , tmp-1000 .0) 

icspec=iconpc(3) 

dl= cpO (icspec ,j c) 

d2= cpl (icspec ,j c) 

d3= cp2(icspec, jc) 

d4= cp3(icspec , j c) 

d5= cp4(icspec , j c) 

icspec=iconpc(l) 

dl=dl+ f lbar(ijk)*cpO(icspec, jc) 
d2=d2+ f lbar(ijk)*cpl(icspec, jc) 
d3=d3+ f lbar(ijk)*cp2(icspec, jc) 
d4=d4+ f lbar(ijk)*cp3(icspec, jc) 
d5=d5+ f lbar(ijk)*cp4(icspec, jc) 
icspec=iconpc (2) 
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dl=dl+ f 2bar(ijk) *cpO(icspec, jc) 
d2=d2+ f2bar(ijk)*cpl(icspec, jc) 
d3=d3+ f 2bar(ijk)*cp2(icspec, jc) 
d4=d4+ f 2bar (ijk) *cp3(icspec, jc) 
d5=d5+ f 2bar (ijk)*cp4(icspec, jc) 
cpbar ( i j k) = ( ( (d5*tmp+d4) *tmp+d3) *tmp+ 

> d2)*tmp+dl 

ktbar (ijk) =mubar ( i j k) *cpbar (ijk) /prandl 

endif 

RETURN 

END 

function find_inverse_molecular.weight.mcs (element , centroid) 
include J cf sparms . i ; 
include ' cf svars . i J 
include * cf sgas . i } 
c include , cfschem < i ; 

integer element 
logical centroid 
c 

include } p3dpar . i ' 
include ; p3dpro . i } 

c 

c 

c PURPOSE: Computes inverse of the molecular weight of a gaseous 
c mixture at a nodal point of the computational grid. 

c 

c 

icspec=iconpc (3) 

f ind_inverse_molecular_weight_mcs= dwm(icspec) 
if(centroid) then 
icspec=iconpc (1) 

find_inverse_molecular .weight _mcs= 

> f ind_inverse_molecular_weight_mcs+f 1 (element )*dwm(icspec) 
icspec=iconpc(2) 

f ind.invers e. molecular .weight _mcs= 

> f ind. inverse. molecular. weight _mcs+f 2 (element ) *dwm(icspec) 
else 

icspec=iconpc (1) 

f ind_inverse.molecular_weight.mcs= 

> f ind. inverse.mol ecular. weight .me s+f lbar (element) *dwm(icspec) 
icspec=iconpc(2) 
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f ind_inverse_molecular_weight_mcs= 

> f ind_inverse_molecular_weight_mcs+f 2bar (element ) *dwm(icspec) 

endif 
C 

return 

end 

function f ind_h_mcs (element ,centroid) 



include 

* cf sparms . i J 


include 

* cf svars . i > 


include 

' cf shref . i * 


include 

} cf sgas . i J 


include 

' cf scpmix . i J 

c 

include 

: ; cfschem.i J 

c 

include 

'BLOCK' 


integer 

element 


logical 

centroid 

c 

include 

; p3dpar . i ; 


include 

; p3dpro . i ; 

c 

c 




c PURPOSE: Computes specific enthalpy at a nodal point of the 
c computational grid. 

c 

c 

if (centroid) then 
tmp = t (element) 
j c=cvmgp(2 , 1 , tmp-600 .0) 
j c=cvmgp(3 , jc , tmp- 1000 .0) 
icspec=iconpc(3) 
hf mm0=dhf orm(icspec) 
dl= cp0(icspec, jc) 
d2= cpl (icspec , jc) 
d3= cp2(icspec, jc) 
d4= cp3(icspec, jc) 
d5= cp4(icspec , jc) 
icspec=iconpc(l) 

hfmm0=hf mm0+ f 1 (element) *dhf orm(icspec) 
dl=dl+ f 1 (element ) *cp0 (icspec , jc) 
d2=d2+ f 1 (element) *cpl (icspec , jc) 
d3=d3+ f 1 (element) *cp2 (icspec , jc) 
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d4=d4+ f 1 (element)*cp3(icspec, jc) 
d5=d5+ f 1 (element) *cp4( icspec , jc) 
icspec=iconpc(2) 

hfmmO=hfmmO+ f 2 (element ) *dhf orm(icspec) 
dl=dl+ f 2 (element) *cpO (icspec, jc) 
d2=d2+ f 2 (element) *cpl (icspec, jc) 
d3=d3+ f 2 (element ) *cp2 (icspec , j c) 
d4=d4+ f 2 (element) *cp3( icspec, jc) 
d5=d5+ f 2 (element ) *cp4 (icspec , j c) 

c calculate enthalpy of mixture 

f ind_h_mcs=( ( ( (0 . 20*d5*tmp+0 . 25*d4) *tmp+d3/3 . 0) *tmp+ 

> 0 . 50*d2) *tmp+dl) *tmp+hfmm0 
else 

tmp = tbar (element) 
j c=cvmgp (2,1, tmp-600 . 0) 
j c=cvmgp(3, jc, tmp-1000 .0) 
icspec=iconpc(3) 
hfmmO=dhform( icspec) 
dl= cpO(icspec, jc) 
d2= cpl (icspec ,j c) 
d3= cp2(icspec, jc) 
d4= cp3(icspec , j c) 
d5= cp4(icspec , j c) 
icspec=iconpc (1) 

hfmm0=hfmm0+ f lbar (element) *dhf orm(icspec) 
dl=dl+ f lbar (element) *cp0 (icspec, jc) 
d2=d2+ f lbar (element )*cpl (icspec, jc) 
d3=d3+ f lbar (element ) *cp2 (icspec , j c) 
d4=d4+ f lbar (element) *cp3( icspec, jc) 
d5=d5+ f lbar (element) *cp4( icspec , jc) 
icspec=iconpc(2) 

hfmm0=hfmm0+ f 2bar (element) *dhf orm(icspec) 
dl=dl+ f 2bar (element) *cp0 (icspec , jc) 
d2=d2+ f 2bar (element ) *cpl (icspec , jc) 
d3=d3+ f 2bar (element ) *cp2 (icspec , jc) 
d4=d4+ f 2bar (element) *cp3 (icspec , jc) 
d5=d5+ f 2bar( element ) *cp4( icspec , jc) 

c calculate enthalpy of mixture 

f ind_h_mcs= (( ( (0 . 20*d5*tmp+0 . 25*d4) *tmp+d3/3 . 0) *tmp+ 

> 0 . 50*d2) *tmp+dl) *tmp+hfmm0 
endif 
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return 

end 
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Appendix X 

An Example Summary of CPU Times Taken By CORSAIR and 

EUPDF 

Table 1 summarizes the cpu times per cycle taken by EUPDF and CORSAIR 
solvers versus the number of processors used on the NASA LeRC LACE clus- 
ter. These computations refer to the case of a confined swirl-stabilized spray 
flame as reported in Ref. 7. The computations were performed on a grid of 
3600 quadrilateral elements and a total of 0.36 million Monte Carlo parti- 
cles ( = 100 particles/cell). Both the PDF and CFD solvers showed excellent 
parallel performance with an increase in the number of processors. For a 
discussion of the parallel performance of these solvers refer to Ref. 7. It 
takes approximately about 1000 to 2000 cycles for the computations to reach 
a converged solution. 


Table 1. Cpu time (sec) per cycle versus number of PEs on LACE Cluster. 


Number of processors 

Solver 

Characteristic 

2 

4 

8 

16 

EUPDF 

1 step/cycle 

2.30 

1.35 

0.75 

0.44 

CORSAIR 

5 steps/cycle 

3.55 

1.90 

1.10 

0.60 
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